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Abstract. In this lecture I give an introduction to the rotational energy extraction of black holes 
1^ . by the electromagnetic Blandford-Znajek process and the generation of relativistic jets. After some 

' basic material on the electrodynamics of black hole magnetospheres, we derive the most important 

results of Blandford and Znajek by making use of Kerr-Schild coordinates, which are regular on 
■ the horizon. In a final part we briefly describe results of recent numerical simulations of accretion 

flows on rotating black holes, the resulting large-scale outflows, and the formation of colUmated 
relativistic jets with high Lorentz factors. 
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S : 1. INTRODUCTION 

For good reasons it is by now generally believed that active galactic nuclei (AGNs) are 
i powered by supermassive black holes (BHs), with typical masses of about lO'^M©. In- 

deed, this seems to be the only way to generate the enormous energy of order lO'^^ erg/s 
OO • in regions not much larger than the solar system. Accretion of matter by black holes, of- 

^ ! fers the most efficient power supply. A few solar masses of gas per year suffice to power 

I the most energetic quasars. 

^ ■ It is known for long that AGNs often expel enormous energies in two oppositely di- 

O ' rected relativistic jets with Lorentz factors 7 ~ 10. More recently, moderately large 

Lorentz factors (7 a few) have been seen in BH x-ray binaries (XRBs), called micro- 
quasars. Relativistic jets have also been observed in gamma-ray bursts (GRBs). How are 
these remarkable jets formed? Which mechanisms are responsible for the concentrated 
5^ \ energy input at their origin? Most probably processes connected to BHs are involved. 

Among these, perhaps the most promising possibility is the energy extraction from a 
black hole via magnetic fields, since fairly strong magnetic fields are likely to be present 
in accretion flows on the central BH. This scenario, which gives rise to a long-range 
coherence, is the main subject of this talk. Basic questions to be answered are: Which 
forces drive the jets? What are the mechanisms that regulate their content, and how are 
they coUimated? 

Penrose first discovered that a substantial fraction of a Kerr black hole mass can 
be converted, at least in principle, into the energy of surrounding matter or radiation. 
However, the Penrose process is inefficient under typical astrophysical conditions. Later 
it was found that external electromagnetic fields can be used to extract rotational energy 
of black holes. Of great influence was a pioneering study by Blandford & Znajek [8] 
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which demonstrated the possibility of an electromagnetically driven wind from a rotating 
black hole, provided the space around the black hole is filled with plasma. In their paper 
they developed a general theory of force-free steady-state axisymmetric magnetospheres 
of black holes, and estimated that the power of the wind could be high enough to explain 
the energetics of radio galaxies and quasars. 

This work triggered the so called membrane model by Thome and collaborators [1] 
(see also [7]). In this model one first splits the elegant 4-dimensional physical laws of 
general relativity (GR) into space and time (3+1 splitting). For a general situation this 
can be done in many ways (reflecting the gauge freedom in GR) since there is no canon- 
ical fibration of spacetime by level surfaces of constant time. (For a stationary BH we 
shall choose foliations which are adapted to the corresponding Killing field.) Relative 
to these the dynamical variables (electromagnetic fields, etc) become quantities on an 
absolute space which evolve as functions of an absolute time, as we are accustomed to 
from non-relativistic physics. We shall see, for example, that the 3+1 splitting brings 
Maxwell's equations into a form which resembles the familiar form of Maxwell's equa- 
tions for moving conductors. We can then use the pictures and our experience from 
ordinary electrodynamics. 

In a second more specific step one replaces the boundary conditions at the horizon 
by physical properties (electric conductivity, etc) of di fictitious membrane. This proce- 
dure is completely adequate as long as one is not interested in fine details very close to 
the horizon. The details of this boundary layer are, however, completely irrelevant for 
astrophysical applications. (The situation is similar to many problems in electrodynam- 
ics, where one replaces the real surface properties of a conductor and other media by 
idealized boundary conditions.) 

This approach turns the drawback of the Boyer-Lindquist coordinates, which become 
singular at the horizon, into an advantage. The membrane model has the merit that it 
allows us to understand astrophysical processes near black holes more easily, because 
things become then closer to the intuition we have gained from other fields of physics, 
for instance from the electrodynamics of moving bodies. However, this membrane 
reformulation of regularity requirements at the horizon - in terms of the singular BL 
coordinates - is clearly artificial. More importantly, it hides the fact that the key role 
in the electrodynamic Blandford-Znajek mechanism is played not by the black hole 
event horizon, but by its ergo-sphere (as in the Penrose process). Later I shall present an 
alternative derivation of the most important results of Blandford and Znajek, in making 
use of Kerr-Schild coordinates. This foliation is not singular at the horizon, whence no 
boundary problems arise. 

Kerr-Schild coordinates not only simplify theoretical analysis, but have also success- 
fully been used in some numerical simulations. In recent years several codes have been 
developed to study the magnetic energy extraction from BHs. In particular, various 
groups have addressed the question how relativistic jets may have formed. In the last 
part of the lecture I will briefly report on what has been done and achieved so far in 
numerical studies of what is now often dubbed the (generic) "Blandford-Znajek mecha- 
nism", although only a specific version was originally proposed. 

We begin with some basic material on black hole electrodynamics. 



ENERGY EXTRACTION FROM BLACK HOLES 



Febraary 2, 2008 



2 



2. SPACE-TIME SPLITTING OF ELECTRODYNAMICS 



Let us perform the 3+1 splitting of the general relativistic Maxwell equations on a 
stationary spacetime {M,^^')g). Most of what follows could easily be generalized to 
spacetimes which admit a foliation by spacelike hypersurfaces (see, e.g., Ref. [2]), but 
this is not needed in what follows. Similar 3+1 decompositions can be carried out for 
the other equations of general relativistic magnetohydrodynamics (GRMHD) . 

Slightly more specifically, we shall assume that globally M is a product R x L, 
such that the natural coordinate ? of R is adapted to the Killing field k, i.e., k = dt. 
We decompose the Killing field into normal and parallel components relative to the 
"absolute space" (E, g), g being the induced metric on Z, 

dt = au + p. (1) 

Here u is the unit normal field and /3 is tangent to E. This is what one calls the 
decomposition into lapse and shift; a is the lapse function and /3 the shift vector field. 
We shall usually work with adapted coordinates (jc'^ ) = , jc') , where {jc'} is a coordinate 
system on L. Let j8 = j8'<9; {di = d/ djd), and consider the basis of 1 -forms 

adt, dx' + P^dt. (2) 

One verifies immediately, that this is dual to the basis {u, df} of vector fields. Since u is 
perpendicular to the tangent vectors di of E, the 4-metric has the form 

(4)g = -a^dt^ + gij {dx' +p'dt) {dxj + pjdt) , (3) 

where gijdx'dx^ is the induced metric g on E. Clearly, a, j8, and g are all time- 
independent quantities on E. 

We introduce several kinds of electric and magnetic fields. Obviously, the spatial 1- 
forms 

i:=-ikF, J^:=ik*F, (4) 
have an intrinsic meaning. If we set 

F = ^ + ^Adt, *F = ^-J^Adt, (5) 

the electric and magnetic 2-forms are spatial. Two further electric and magnetic 

1 -forms are defined by 

<f = *^, J^ = *^, (6) 

where ★ denotes the spatial Hodge dual on (E,g). The corresponding vector fields are 

denoted by E and B. These are measured by observers moving with 4-velocity w, so- 
called FIDOs, for fiducial observers. One easily finds the following algebraic relations: 

i = a^-ip^, je' = aM' + ip^. (7) 

We give also the coordinate components of the various fields: 

ii = -Fti, M'i = ^Fn = ^nijkF^\ ^ij = Fij, = ^r]'j^Fjk, ^ij = *Fij, E' = aF'^ 

(8) 
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where riij^. denotes the Levi-Civita tensor on (E, g). 

In terms of the current 4-vectory^ the electric charge density is p^i = ocJ' . Beside the 
spatial current 1-form j — Jj^dx'^ we use also its spatial Hodge dual ^ —-><]■ Using also 
P = *PeZ, the Maxwell equations can be written in 3+1 form as (d denotes the exterior 
differential onL): 

d^ = 0, dt^ + Ai = Q, 
d^ = 47tp, -dt^ + dJ^ = 47t/, J ■.= a/ -ipp. (9) 

These can be translate into a vector analytic form. (For this as well as detailed deriva- 
tions, see [7].) All this looks like Maxwell's equations for moving conductors. 

Integral Formulas 

As is well-known from ordinary electrodynamics, it is often useful to write the basic 
laws in integral forms. Consider, for instance, the induction law in (9). If we integrate 
this over a surface area which is at rest relative to the absolute space, we obtain with 
Stokes' theorem := ds^), using also (7), 




The left hand side is the electromotive force (EMF) along ^. The last term is similar to 
the additional term one encounters in Faraday's induction law for moving conductors. 
It is an expression of the coupling of ^ to the gravitomagnetic field and plays a crucial 
role in much that follows. This term contributes also for a stationary situation, for which 
(10) reduces to 

EMF(<^) ^ I aS ^ I (11) 

J ^ J ^ 

3. BLACK HOLE IN A HOMOGENEOUS MAGNETIC FIELD 



As an instructive example and a useful tool we discuss now an exact solution of 
Maxwell's equations in the Kerr metric, which becomes asymptotically a homogeneous 
magnetic field. This solution can be found in a strikingly simple manner [3]. 

For any Killing field K and its 1-form one has the following identity 5dK^ = 
2R{K), where 5 denotes the co-differential and R{K) is the 1-form with components 
RnyK^. In components this is equivalent to 

^^"a = —R^aK". (12) 

This form can be obtained by contracting the indices a and p in the following general 
equation for a vector field ^a;pn — ^a;np — ^x^apn by using the consequence 
= of the Killing equation K<j-p +Kp-(j = 0. For a vacuum spacetime we thus have 

ddK^ = for any Killing field. Hence, the vacuum Maxwell equations are satisfied 

if F is a constant linear combination of the differential of Killing fields (their duals, 
to be precise). For the Kerr metric, as for any axially symmetric stationary spacetime. 
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we have two Killing fields k and m. Both in Boyer-Lindquist (BL) and in Kerr-Schild 
(KS) foliations the Killing fields are coordinate derivatives: k — dt and m = d(p (k^ and 
denote the corresponding 1 -forms). The Komar formulae provide convenient expressions 
for the total mass M and the total angular momentum / of the Kerr BH: 



(for G = 1). 
We try the ansatz 



M=-^ [ *dk\ J=T^ I *dm' (13) 
StuJoo ' 1671:' ^ ^ 



F=l-Bo {dm^ + la dk^ ) {Bq = const) , (14) 



2 

and choose a such that the total electric charge 



An 

vanishes. The Komar formulae (13) tell us that 

Q = -^Bo(l67i:J-2a-S7i:M), (16) 

and this vanishes if a = J/M (which is the standard meaning of the symbol a in the Ken- 
solution). 

Clearly, F is stationary and axisymmetric: 

LkF=LrnF = 0, (17) 

because (dropping b from now on) L/^dk = dLj^k = 0, {L]^k = [k,k] = 0), etc. 
The solution (14) can be expressed in terms of a potential: F = dA, with 

A^^Bo{m + 2ak) = ^Bq {g^^ + lag^t) dx^ , (18) 

where the last expression holds in Boyer-Lindquist (BL) as well as in Kerr-Schild (KS) 
coordinates. Asymptotically, this describes a magnetic field in the z-direction whose 
magnitude is Bq. 

It is straightforward to work out explicit expressions for E and B. The electric field 
has a quadrupole-like structure and is poloidal. It is proportional to a, and thus due to the 
gravitomagnetic component of the Kerr solution. Its emergence is of great astrophysical 

significance. 

It is of interest to work out the magnetic flux through the equator of the BH, i.e., 

^= / [ , ^ = 271^^ equator- (19) 

jupper h. Jequator ^ 

One finds 

<i> = 47tBoM{rH -M) = 47tBoM\/M^-a^. (20) 
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Generically one has, as expected, <I> ^ tzv^Bq. Note, however, that (20) vanishes for an 
extremal BH (a = M). In other words, the flux is completely expelled from the black 
hole, like in the Meissner-Ochsenfeld effect in superconductivity. In order to study the 
structure the Wald solution close to the black hole one has to express it in terms of Kerr- 
Schild coordinates (see Sect. 6). The magnetic field lines in these coordinates are shown 
for a = M in Fig. 1, taken from [4]. 




-2 2 



X 

FIGURE 1. Black hole "Meissner effect" (from [4]). The magnetic field lines of Wald's vacuum 
solution in KS coordinates for a = 1 are shown in the x = r sin tJ, z = rcos -& plane. The partial circles 
(thick lines) bound the ergosphere. 



4. AXISYMMETRIC STATIONARY FIELDS 



In this section we discuss some consequences of Maxwell's equations for a stationary 
axisymmetric magnetoshere outside a black hole. 



4.1. Potential Representation 

For an axisymmetric field we have L^^ ^ = ^ dig^ = 0, whence i^^ SS = 
—d^ 12%. From this we conclude that for the BL foliation can be expressed in terms 
of two potentials ^ and /, 

^=-L(i^Ad(p+ — *d(p , (21) 

poloidal part toroidal part 

both of which can be taken to be independent of (p. (The physical meaning if / will 
be discussed further below.) ^ is the magnetic flux function (see Fig. 2), because the 
poloidal flux inside a tube = const} is 

j j d{^'d(p) = ^j^'d(p = ^', ^(0)=0. (22) 
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^ is constant along magnetic field lines, as should be clear from Fig. 2. It is easy to 




FIGURE 2. Axisymmetric magnetic field. The total flux inside the magnetic surface defines the flux 
function ^. 

show this also formally. 

The electric field E has no toroidal component for an axisymmetric stationary sit- 
uation. This is an immediate consequence of the induction law: Applying (10) for a 
stationary and axisymmetric configuration to the closed integral curve ^ of the Killing 
field d(p, we obtain 

^•«e^^= / d^ = 0^£^°^ = 0. (23) 

A similar application of Ampere's law in integral form gives 

I = ATt f a J =AkI, (24) 

where / is the total upward current through a surface bounded by ^ . This shows that 
the potential / in (21) is the upward current. 

So far we did not make any model assumptions about the physics of the magne- 
tosphere. But now we assume that the electromagnetic field is degenerate, i.e., that 
E ■ B = 0, which is equivalent to the invariant statement F AF = 0. (This is satisfied 
for ideal GRMHD and for force-free fields; see Sects. 5 and 6.) Because S" is poloidal, 
we can then represent the electric field as follows 

^ = k^^ {E = -vfxB), (25) 
where Vf is toroidal. Let us set 

1 . 

Vf ='■ —{^F-co)COe(p. (26) 

For the interpretation of Q.f note the following: For an observer, rotating with angular 
velocity Q., the 4- velocity is u = u' {dt + Q. dq,) . On the other hand, u = Y{eo + v), where 

V is the 3-velocity relative to a FIDO. Using also dt = aeo + P we get Q.d(p = av — P or 

v = — {Q.-co)d(p = —{Q.-co)a)e(p. (27) 
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This has the same form as (26). Since the transformed electric field E' = Yf{E + vf x B) 
vanishes, we may regard Q.f as the angular velocity of the magnetic field lines. (The 
transformed frame can be defined as the "local rest frame" of the magnetic field lines.) 
Eqs. (25) and (26) imply 

aS = d^, ^ = -d^, (28) 

thus E is perpendicular to the surfaces = const}. Taking the exterior derivative, and 
using induction law, we get A£If A = =^ = ^f(^)- So the electromagnetic 
field is determined in terms of ^, / and the "flow (field line) constant" ^p{^) . 



4.2. EMF outside a rotating Black Hole 



In Fig. 3 we consider a stationary rotating BH in an external magnetic field (like in 
§3). The integral in (1 1) along the field lines gives no contribution and far away /3 drops 
rapidly (~ r^^). Thus, there remains only the contribution from the horizon (^//) of the 
path in Fig. 3: 

^H = -^Hd^, (29) 



where ^2//is the angular velocity of the horizon (only the normal component 5^ con- 
tributes). I recaU that Q.h = a{2MrHy\ rH=M + VM^-a^. 




FIGURE 3. Arrangement for eq. (29). 

Let us work this out for the special case of an axisymmetric field. For the closed path 
<^ in Fig. 3 the EMF is by (29) 



EMF = AV 



H 



lit 



d^, 



I.e. 

Ay = — A^. 

In 

This result is independent of the physics outside the black hole. 
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Let us integrate it from the pole to some point north of the equator. For the exact 
vacuum solution in §3 we know the result for the EMF, if we integrate up to the equator: 
From (30) and (20) we get EMF = ISIhBqM (rn-M) or {Qh = a/lMrn) 



EMF = aBo ^ (rH = M+\/M^ -aA . (31) 



Note that the "Meissner effect" for black holes implies that the EMF vanishes for ex- 
tremal black holes. We shall, however, see that this property of the vacuum solution is 
astrophysically not relevant, because a plasma-filled magnetosphere changes the struc- 
ture of the magnetic field close to the horizon dramatically. Even for a maximally rotat- 
ing Kerr black hole the magnetic field is pulled inside the event horizon. For a detailed 
discussion we refer to [4]. 

2 

For a general situation we have roughly E = ^ B^Tir^, a)^~<ft)^>~^. 
The total EMF, V = LAV, is thus 

1 la n I / a \ 

Vr^—^H'Vr^———B^nrjjC^-{ — )MB^_ (32) 

(compare this with (31)). Numerically we find 

y ~ (10^0 Volt) f^)^^^ (33) 

^ ^ \m) 109 M© 104 g 

For reasonable astrophysical parameters we obtain magnetospheric voltages V ~ 10^^ 
Volts. This voltage is comparable to the highest cosmic ray energies that have been 
detected. 

Note, however, that for a realistic astrophysical situation there is plasma outside the 
BH and it is, therefore, at this stage not clear how the horizon voltage (33) is used in 
accelerating particles to very high energies. This crucial issue is addressed in the final 
section. 

Let us estimate at this point the characteristic magnetic field strength than can be 
expected outside a supermassive BH. A characteristic measure is the field strength 
Be, for which the energy density B\/%k is equal to the radiation energy density ue 
corresponding to the Eddington luminosity. One finds (Afe g = M/z/lO^M©) 

5£ = 1 .2 X 10^ M~ Gauss. (34) 

For a BH with mass ~ 10^ M© inside an accretion disk acting as a dynamo, a character- 
istic field of about 1 Tesla (10** Gauss) is thus quite reasonable. 



5. BASIC EQUATIONS OF GENERAL RELATIVISTIC IDEAL 
MAGNETOHYDRODYNAMICS 

The relativistic fluid is described by its rest-mass density, po, the energy-mass density, 
p, the 4-velocity, U^, and the isotropic pressure, p, assumed to be given by an ideal gas 
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equation of state 

P=(r-l)e, (35) 

where e = p — po is the internal energy density, and F is the adiabatic index. 

The basic equations of GRMHD are easy to write down. First, we have the baryon 
conservation 

V^(poC/^)=0. (36) 
For a magnetized plasma the equations of motion are 

Vyr^'' = 0, (37) 

where the energy-stress tensor T^^ is the sum of the matter (M) and the electromagnetic 
(EM) parts: 

= {p+p)U^U'+pg^\ (38) 
Te'm = l^iF^xF'''-\g^''F,pF-P). (39) 

In addition we have Maxwell's equations 

dF = 0, VyF^^^AnJ^. (40) 

We adopt the ideal MHD approximation 

iuF = 0, (41) 

which implies that the electric field vanishes in the rest frame of the fluid (infinite 
conductivity). Then the inhomogeneous Maxwell equation provides the current 4-vector 
J^, but is otherwise not used in what follows. 

As a consequence of (40) and (41) we obtain (using the Cartan identity Lu = iu o 
d + doiu): LuF = 0, i.e., that F is invariant under the plasma flow, implying flux 
conservation. The basic equations imply that 

*F = hAU, (42) 

where h = i[j*F is the magnetic induction in the rest frame of the fluid (seen by a comov- 
ing observer). Note that iuh = 0. Furthermore, one can show that the electromagnetic 
part of the energy-momentum tensor may be written in the form 



(43) 



with \\hf ^ hah". 

Let us also work out the 3+1 decomposition of the ideal MHD condition (41). Using 
U = Y{eo + v), 7 = ( 1 - v^) ~ and the coordinate velocity V = U'/ C/^ = av' - j8 we 
obtain from (5) 

^ = i^^ or i = i^^ (^i^i^O). (44) 
Therefore, the induction equation becomes 

dt^ + di^^ = or dt^ + L^^ = 0. (45) 
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6. THE BLANDFORD-ZNAJEK PROCESS 



In this section we derive the main results obtained by Blandford and Znajek [8], by 
making use of the Kerr-Schild coordinates. 



6.1. Steady-state force-free magnetospheres 

BZ studied axisymmetric, force-free magnetized plasmas outside Kerr black holes. 
The following presentation is strongly influenced by the recent paper [9]. 
A plasma is said to he- force-free if ijF — 0, i.e. 

F^vr = (46) 

(no electric field in the rest system of the current). This condition follows from ideal 
GRMHD when the inertia of the plasma is ignored. Formally, it is obtained by letting the 
specific enthalpy (p +p)/n in the ideal GRMHD equations go to zero. Since Maxwell's 
equations imply that "^vT^m ~ —F^^P, the energy-stress tensor of the electromagnetic 
field is separately conserved. The force-free condition (46) implies the constraint F A 
F = 0, i.e., the algebraic condition *Ffxy,F^'^ = 0. In the literature it has often been 
asserted that magnetospheres of black holes should in large parts be nearly force-free 
This is not born out in recent simulations, except in the polar region (see Sect. 7). In 3-1-1 
decomposition (46) becomes 

peE + jxB = 0. (47) 
Therefore, the vector fields E and B are perpendicular. 



6.2. Energy flux at infinity 

Of particular interest is the energy flux at infinity. Since dt is a Killing field, T^^ are 
the components of a conserved 4-vector field: 



1 



-8 



For stationary fields ^i{^/^T/) — 0. The electromagnetic power at infinity^ is 

Pem= r V^FslTtd^, (48) 

^0 



In the classical BZ process one ignores the extraction of rotational energy of the black hole by accreting 
material interacting with the electromagnetic field. Below T^^^ always denotes the electromagnetic energy- 
momentum tensor. 
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where Fe = —T^. The integral can be taken for any fixed r outside the horizon, in 
particular at the horizon if we use Kerr-Schild coordinates. For Fe we get 

Fe = -—FtiF'' = - —Ft^F'^ (49) 

(Ft(p —0). It is straightforward to check that Fe is invariant under the coordinate trans- 
formation [BL]— >^ [KS]. Using previous results on readily finds in BL coordinates 

Fe = -^aF-^^^B'B^, (50) 

where 5' :^B'/a = 

Now we transform the result (50) to KS coordinates. It is straightforward to show that 
B^ remain invariant, while the remaining component transforms as 

B^iBL] =B^[KS]-^-^^B'^[KS], (51) 



where A := — 2Mr + a^. This implies that on the horizon (A = 0) 

1 

2n 

Therefore, 



Fe\h = :^ {B'f^FrHi^H - ^f) sin^ ^. (52) 



Pem = j di»pjjSm'&{B'')^ QprHi^H-^F) sin^^ (53) 
or in terms of the normalize normal component B± 

Pem = \ fd^QlFiO^H - 0-f)Ph(0^b\. (54) 

2 Jo 

The important result (54) shows that Pem becomes maximal foiQ.F ^ k^H- Then 

P^M^l^lf pH&^Bld^. (55) 
After angular integration one finds, if B\ is replaced by an average value 

where the function / is not far from 1 [10] . Numerically, 

P'SS! = 1.7 X 10^6 ^ f ^ ■ . ^V/. (57) 

For the angular momentum flux Fi one finds Fe = ^fFl- Using this one can show 
that up to 9% of the initial mass can, in principle, be extracted [10]. 

Whether the crucial condition Q.e ~ is approximately satisfied in realistic as- 
trophysical scenarios is a difficult problem for model builders. This brings me to recent 
numerical work by several groups. 
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7. GENERAL RELATIVISTIC MHD SIMULATIONS 



In recent years several groups have developed codes for ideal GRMHD, and applied 
them in particular for studies of the generic BZ mechanism. These numerical studies 
show how accretion dynamics self -consistently create large scale magnetic fields and 
explore the resulting outflows. 

7.1. Numerical methods 

There are many ways to write the basic equations in a form suitable for numerical 
integration. A pioneering code, based on the 3+1 splitting relative to the FIDO tetrad 
described in Sect. 2, has been developed by Koide and collaborators (for a detailed de- 
scription, see [11]). The system to be integrated consists of eight evolutionary equations 
for two scalar and two vector quantities (primary code variables), which contain in ad- 
dition five 'primitive' variables that can be determined from the former by solving two 
coupled polynomial equations P{x,y) = 0, Q{x,y) = 0. This process, like similar itera- 
tive root-finding processes in other schemes, is time consuming. 

A code that is able to perform long-term (several thousand of M in time) GRMHD 
simulations, has been developed by Villiers & Hawley [12]. It evolves different auxiliary 
variables from which the primitive variables are easily recovered. In contrast to some 
alternative schemes [13], [14], it is not fully conservative. Since BL coordinates are used, 
the inner boundary condition must lie outside the horizon. KS coordinates are used in 
the axisymmetric code HARM [13]. In a modified version of this [19] the inversion from 
"conserved" quantities to "primitive" variables has been improved (see Appendix A). 

A specific common problem is that the magnetic field is advanced in time by an 
antisymmetric differential operator, and not by a differential operator of divergence 
form. In addition one has to guarantee that the numerical scheme preserves the constraint 
equation = to rounding error. One way to handle this problem is the so-called 
constraint transport (CT) method, where the induction equation is discretized such that 
the solenoidal constraint is built in [15]. 

No code is perfect, and it is therefore important that the outcomes of different ap- 
proaches are compared. 

7.2. Qualitative results 

I first summarize some of the results presented in [16] on simulations of accretion 
flows on rotating black holes, and the properties of the resulting unbound outflows^. 

The qualitative late-time structure is illustrated in Fig. 4. Along the equator there 
is a wedge-shaped Keplerian accretion disk and a net accretion flow is produced by 
MHD turbulence whose origin is a magneto-rotational instability, which leads to a 



This is based on a sequence of earlier papers on the the subject, cited in [16] 
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fast amplification of the magnetic field (see Appendix B). A short distance outside the 
marginally stable circular stable orbit, the equatorial pressure and density reach a local 
maximum in a 'inner torus' region. Inside this local pressure maximum, the density and 
pressure drop as the flow accelerates toward the black hole in the 'plunging region'. 
Above and below the disk is a 'corona' of hot magnetized plasma with a magnetic field 
whose typical strength is near equipartition. Along the spin axis of the BH there is a 
'centrifugal funnel' that is largely empty of matter, but filled with magnetic field and an 
outward Poynting flux. Between the evacuated funnel and the corona there is a region 
of unbounded mass flux, referred to as the 'funnel wall jet'. It is this part that is the 
center of attention in [16]. Of particular interest is the strength of the jet compared to 
the amount of accretion on the BH. The quantification of this is somewhat ambiguous, 
because matter and electromagnetic fluxes are usually not constant in time or radius, and 
there is exchange between the two components. 




jet 



FIGURE 4. Main dynamical features in accretion disk simulations. Fig. 1 of [16]. 

The numerical results show that for rapidly-rotating BH's the jet plays a significant 
role in the energy budget. When | a/M | > 0.9, the total jet efficiency is generally several 
tens of percent, with the matter portion somewhat larger than the electromagnetic part. 
The ratio of angular momentum in the unbound outflow to that deposited in the BH is 
also a very strong function of BH spin. In the non-rotating case this is only a few tenths 
of a percent, but rises to ~ 25% when a/M = 0.9. Vox a/M = 0.99, the highest spin case 
reported in [16], the rate at which the BH spins down due to electromagnetic torques 
nearly matches the rate at which it acquires angular momentum by accreting matter. 
Indeed, the angular momentum expelled in the outflow is an order of magnitude larger 
than the net amount captured. In all cases, the electromagnetic portion of the angular 
momentum carried away is comparable to the matter portion carried in the funnel wall 
jet. 

It has to be emphasized that in the initial state of the simulations there is no mag- 
netic field in the region that eventually becomes the outflow region. The large scale 
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magnetic field witliin tlie funnel and the Poynting flux jet rapidly develop as a result 
of magnetically controlled accretion. The formation can be described as follows. When 
disk material reaches the horizon, strong magnetic pressure gradients are built up which 
drive the plasma upward. In turn, this motion drains off the field lines, forming a "mag- 
netic tower". This was previously seen in non-relativistic simulations using a pseudo- 
Newtonian potential [18]. It is interesting that the evacuated funnel is the only region 
that is force-free^ (in contrast to earlier expectations). 

The matter-dominated outflow moves at a modest velocity (v/c ~ 0.3) along the 
centrifugal barrier surrounding the evacuated funnel. The funnel wall jet turns out to 
be accelerated and collimated by magnetic and gas pressure forces in the inner torus 
and the surrounding corona. The magnetic field is spun by the rotating Kerr spacetime, 
hence the energy of the Poynting flux jet comes from the BH rotation. Below we will 
address the question whether magnetic forces might provide additional acceleration and 
coUimation on far larger scales than are modeled in [16]. 

Hawley and Krolik conclude from their numerical studies that while the proximate 
energy source is the BH's rotation, accretion replenishes both the BH's mass and angular 
momentum. A substantial decrease of the rotational energy of the BH does not appear 
to be a generic phenomenon. This is in contrast to the classic BZ model. For a detailed 
discussion, I refer to the original paper. 

Similar results were found before by McKinney and Gammie in their axisymmetric 
simulations [9]. Among other aspects these also showed that as the hole loses energy and 
angular momentum, its total mass and angular momentum are replenished by accretion. 
In [20] McKinney followed the evolution of the jet till t ^ XQ^GM/c^ out to r 
lO'^GM/c^, and found that by then the jet has become superfast magnetosonic and 
moves at a Lorentz factor of about 10. This may, however, only be a small fraction 
of the Lorentz factor at much larger distances. Indeed, if a large fraction of the magnetic 
and thermal energy would go into particle kinetic energy, one would estimate from the 
simulation that the terminal Lorentz factor may reach almost 10^. This estimate is based 
on the following fact. For a stationary axisymmetric flow the energy (momentum) flux 
per unit rest-mass flux is conserved along flux surfaces and can thus be determined from 
local flow quantities. The relevant formula for this quantity is, using previous notation, 

PqU^ Po 47t PqU^ 

Since the simulation in [20] stopped long before the end of the acceleration period, 
Komissarov et al. have studied numerically the further evolution of the relativistic jet 
[21]. Key issues of this investigation are: (i) Is the magnetic driving mechanism able 
to accelerate outflows to high Lorentz factors with high efficiency over astrophysically 
extended scales? (ii) Can these flows be collimated by purely magnetic stresses? Since 
most of the acceleration takes place far away from the black hole, the simulations are 
carried out in the framework of special-relativistic ideal MHD. 



' The criterion for the vahdity of this property is that {\\h\\^/Sn)/{p +p)^ 1. 
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The authors investigate models of the following kind. The free boundary with an am- 
bient medium is replaced by solid rigid walls on which appropriate boundary conditions 
are imposed. This simplification enables higher numerical accuracy. At the inlet bound- 
ary, the injected poloidal current distribution is prescribed through a rotational profile. 
The initial configuration corresponds to a non-rotating purely poloidal magnetic field 
with nearly constant magnetic pressure across the funnel. Moreover, the outflows are ini- 
tially Poynting flux-dominated. The authors find that these approach a steady state with 
a spatially extended acceleration region. Furthermore, the acceleration process turns out 
to be very efficient; almost 80% of the Poynting flux is converted into kinetic energy. 
The results also show efficient self-coUimation. In contrast to [20], no instabilities or 
shocks are found in the simulation. 

We have already remarked earlier that even for very rapidly rotating black holes the 
Blandford-Znajek process can drive magnetized jets. Fig. 5 from [4] shows clearly that 
there is no "Meissner effect" at work. This is due to the fact that within the ergosphere 
the plasma unavoidably corotates with the black hole. 




X 



FIGURE 5. Magnetic field lines and logarithm of rest-mass density (colored) quasi-steady accretion 
disk simulations for an almost extreme Kerr black hole {a/M = 0.999) [Fig. 3 of [4]]. 

One of the main shortcomings of existing simulations is that the radiation field is 
completely neglected. It is likely that radiative and other high energy processes play a 
significant role in the flow dynamics through radiation force on the outflowing plasma. 
Moreover, at some point physical resistivity will also have to be included, and the single- 
fluid approximation will not always hold everywhere. Numerical calculations will be 
performed in 3D. 
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A. NUMERICAL SCHEMES FOR GRMHD 



Most GRMHD codes have adopted a conservative scheme, which means that the inte- 
grated evolutionary equations are of the general form 

dtU{P) + diF\P)^S{P). (59) 

Here, U_ are vector- valued "conserved" variables, P is a vector of "primitive" variables 
(rest-mass density, internal energy density, velocity components and magnetic field 
components). The fluxes F', U_ and the "source vector" S depend on P. Conservative 
numerical schemes advance U_, then calculate P{U_) once or twice per time step. 

While in praxis the map P i-^ U_{P) is analytically known, the inverse map U_ t-^ P{U) 
is not available in closed form and must be computed numerically. How this is performed 
is at the heart of a conservative scheme, since the operation must be accurate, fast and 
robust. 

Below we describe these points in more detail, following [19]. We use the notation 
and basic equations introduced in Sect. 5. 

Equations (36), (37) and (45) can readily be written in the form (59) in terms of the 
eight conserved variables 

U: D-.^YPo, Q^:=-UyT\, B, (60) 
where Y = —u^U^ = {1 —v^)"^^^. As eight primitive variables we use 

P : Po, e, B, v\ (61) 

(The vector B is common to U_ and P.) 

In [ 1 9] it is shown that W := pohf^, where h:= l+e+ p/ po, and can be determined 
from t/ by solving two polynomial equations. After that one can easily recover all 
primitive variables. 

B. MAGNETO-ROTATIONAL INSTABILITY 

This important instability, whose crucial astrophysical implications have been under- 
stood astonishingly late, shows up already in linearized theory. For a detailed pedagogi- 
cal discussion we refer to [17]. 

In a systematic treatment one linearizes the basic equations of MHD for small fluc- 
tuations of a disk system, consisting of a central point mass and a differentially rotat- 
ing magnetized disk. An analysis of the resulting rather complicated dispersion relation 
leads to the following results. In the non-rotating limit (homogeneous unperturbed situ- 
ation) one finds the familiar Alfven waves and two other MHD modes. The fast one is 
often referred to as magnetosonic wave, and represents magnetic and thermal pressure 
in concert. In the slow mode magnetic tension and gas compression act in opposition. 
For weak fields it becomes degenerate with the Alfven mode, while for strong fields it 
becomes an ordinary sound wave, channeled along the field lines. 
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The effect of Kepler rotation on these three modes is very interesting: At a critical 
rotation frequency the slow MHD mode becomes unstable. 

There is a simple way to understand this far reaching instability. Consider an axisym- 
metric gas disk in the presence of a vertical magnetic field, that has no effect on the disk 
equilibrium. Assume that a fluid element is displaced from its circular orbit by oc e'^^ 
(z = vertical direction). Using the induction law and simple mechanics one finds the 
these incompressible planar displacements satisfy the same equation as the separation 
of two orbiting mass points, connected by a spring (with a spring constant related to the 
Alfven velocity). It is quite obvious that this system is unstable; the separation of the two 
mass points rapidly increases. This is the essence of the weak-field magneto-rotational 
instability. 

The magneto-rotational instability plays a fundamental role in disk accretion, because 
it leads to disk-turbulence and corresponding stresses. 
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